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Imaginary- and real-time versions of an equation for the condensate density are presented which 
describe dynamics and decay of any spherical Bose-Einstein condensate (BEC) within the mean 
field approach. We obtain quantized energies of collective, finite amplitude radial oscillations and 
exact numerical instanton solutions which describe quantum tunneling from both the metastable 
and radially excited states of the BEC of 7 Li atoms. The mass parameter for the radial motion is 
found different from the gaussian value assumed hitherto, but the effect of this difference on decay 
exponents is small. The collective breathing states form slightly compressed harmonic spectrum, 
n=4 state lying lower than the second Bogoliubov (small amplitude) mode. The decay of these 
states, if excited, may simulate a shorter than true lifetime of the metastable state. By scaling 
arguments, results extend to other attractive BECs. 
PACS number(s): 03.75.Fi, 05.30.Jp, 05.45.-a, 47.20.Ky 

I. INTRODUCTION 

In a series of experiments [fj]-|}, a nonuniform Bose-Einstein condensate (BEC) of 7 Li atoms was formed and proved 
metastable for atom numbers N < N c with N c w 1300, in agreement with theoretical predictions |^-@]. Due to 
attraction between atoms, such condensate is bound to collapse when N > N c , but it may collapse also for N < N c 
via quantum or thermal tunneling. Within the mean-field appraoch, BEC is described by one wave function. Our 
aim in this work was to find the exact mean-field description of the quantum tunneling in this simplest conceivable 
many-body system. Our additional motivation is the recently demonstarted ability to control the interaction of 85 Rb 
atoms in BEC 0, which opens a new perspective to systematic experimental checks on quantum tunneling. 

Up to now, studies of the quantum tunneling of BEC relied on assuming gaussian wave functions ||, at least in 
assigning the mass parameter H . Strictly speaking, once the mean field equation is specified, there is no place for such 
an assumption: This equation, taking a form of the nonlinear field equation, by itself determines the dynamics. To 
find solutions which correspond to quantum decay we use the method of instantons, i.e. fields evolving in imaginary 
time ||, carried over to mean- field theories of many-body systems 0,0. It gives the decay rate T — Ae~ s , with 
the exponent S being the action for the optimal mean-field instanton, called bounce. We find the exact instantons for 
spherical BEC by first transforming the original imaginary-time mean-field equation to an equation for the condensate 
density, and then solving it numerically. Having done that wa can check the gaussian ansatz. 

The real-time version of the obtained equation encompasses finite amplitude collective radial oscillations of BEC. 
Applying quantization rule we find energies of radial eigenmodes. By using imaginary-time dynamics we also find 
periodic instantons determining decay exponents of these breathing modes. In this way, the collective dynamics of an 
attractive spherical condensate close to instability is obtained from the mean field equation. 

We assume that the dynamics of BEC is governed by the time-dependent Gross-Pitaevskii equation (CPE) [ jjjj 
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where Vtrap is the static trap potential and g = Airh 2 a/m, with a the s-wave scattering length and m the atomic mass. 
The wave function is normalized as J d 3 r | ip(r,t) \ 2 — N, with N the total number of atoms in the condensate. 
In the present work we consider a harmonic, spherically symmetric trap Vtrap — \rnuj^r 2 . This suggests choosing 

do = \Jfif mwo as the unit of length, 1/ujq as the unit of time and Ticuq as the unit of energy. We also change the 
normalization of the wave function to unity, Air J | tp(r,t) \ 2 r 2 dr = 1. To simplify equations we work with the 
function (j)(r,t) = rtp(r,t), for which the GPE reads 



where K = 4TrNa/d , and the generalised density p(r,t) — r 2 | i/j(r,t) \ 2 . For a stationary state, 4>(r,t) = 
4>(r)exp(—iet), e being the single-particle energy or chemical potential. For each N < N c there are two station- 
ary states, one metastable and another unstable, at the top of the energy barrier M (cf Fig. 1). 
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II. QUANTUM TUNNELING AND EQUATIONS FOR CONDENSATE DENSITY 



Quantum tunneling of BEC is described by a specific solution to the equation 

dcj) 1 d 2 4> ,1 2 , T _ p . , , , . 
or 2 or^ 2 r z 

obtained from (g) by a transition to imaginary time, £ — * — it [ p~0| , p~3 1 - Now, the density p(r, r) = </>(r, — r)*</>(r, t) 
as <f>(r,t)* is replaced by (j)(r,—T)* upon i — > —it. This makes Eq.(B) nonlocal in time. Bounce has to satisfy the 
boundary conditions of 1) periodicity, <f>(r,T p /2) = cf>(r,—T p /2) — (f>o[r) — r^j (r), with ipo( r ) the amplitude of the 
metastable state, and 2) barrier penetration, i.e. <fi(r, r = 0) has to be some state of BEC at the other side of the 
potential barrier. Eq. (|3|) preserves both the normalization Air J drp(r, r) = 1 and the energy 



1 d^j-r)* dcb(r) 1 2 Kp 2 

2 dr dr 2 PT 2 r 2 ' 



4tt / dr{- yv a y ' }, (4) 



with = TVS the energy of the metastable state. For a bounce starting from the metastable state the period r p 
extends to infinity ]lO||l^Jll| ] . The decay exponent reads p0| 

5 = 47riV / ^ dr f dr^(-r)*|^(r). (5) 
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Since Eq.(^) is real and the boundary value (f)o{r) may be taken real, we assume real <f>(r,T) in the following. 

Now the point is to transform the nonlocal in time instanton equation (|^) into an evolution equation for the 
condensate density. A transformation of the real-time GPE (|^) to a fluid-dynamic form provides an analogy, but is 
conceptually simpler. 

The bounce equation (^), with the boundary conditions specified above, splits into two time- local equations for the 
time-even density p and the time-odd current j(r, r) = — 1/2(0(— r)9 r 0(r) — 0(r)<9 r 0(— r)) 

^ + *=°' (6) 



where the kinetic energy density = d r ^(—T)d r 4>(r) = [l/4(9 rj o) 2 — j 2 }/ 'p. When p is non-negative (which is very 
plausible in the present case, but not granted in general as p(r, r) = <f>(r,T)<j>(r,—T) not <fi(r, t) 2 ), one can define 
a regular, time-odd function % = — i(ln0(r) — ln0(— r)) which allows decomposition = ^J~pe~ x . From this, the 
fluid-dynamic representation of (H) follows, with the velocity field d\jdr = j / p. However, even for arbitrary p, one 
can eliminate j from Eqs.(^,^), which is more convenient. Introducing /(r, r) = j£ p(r' ,r)dr' , so that j = —df/dr 
and /? = df/dr, we automatically fulfil (H), and (0) transforms to the equation for the primitive of the bounce density, 
basic for the imaginary-time dynamics of spherical BEC-s: 

4 |/ + a (i (S^ ) + a lr + K a ( l )] = a 

or 

Notice, that the finite amplitude oscillations of BEC around the metastable minimum are described by the real-time 
version of Eq.(||), in which d 2 f and (d T f) 2 are replaced by — dff and —(d t f) 2 , respectively. 

An alternative, global approach to quantum tunneling is to look for a minimum of action (||) under the condition 

of constant energy ([|) and norm. Indeed, for a regular Xi i- e - positive p, S = 4irN J^ P J^ 2 ^ T Jq° drj 2 /p. Using Eq.(|4|) 

and introducing an observable Q uniquely labelling states along the barrier, explicitly Q(t) = (r 2 ) /N — Air J °° drpr 2 , 
we obtain the following functional 

Q(r P /2) 



S[f] = 2N / dQy/2B(Q)(V(Q)-£) (9) 

Jq(o) 

which is minimized by the primitive of the bounce density (note, that Q(0) < Q(t p /2) for BEC). The potential energy 
V(Q) = V[p(Q)l 
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and the effective mass parameter B(Q) = B[f(Q)] 

poo (^-) 2 

B(Q) = 4tt / dr-^— (11) 
Jo P 

are both the functionals of /. Note that Eq.(||) is invariant with respect to a change of the controlling variable, 
as for any other such variable q, B(q) = B(Q)(dQ/dq) 2 . The energy conservation (Q) implies that for bounce 
Q 2 = 2(V(Q) - £)/B{Q), with Q = dQ/dr, and therefore Q = &[{V{Q) - £)/B{Q)\. 



III. RESULTS AND DISCUSSION 

We have solved Eq.(|J) using the variable Q rather than r. An initial sequence of densities p s (r,Qi), with 30-50 Qi 
points covering the barrier region, was constructed by minimizing V[p] ([l0|) under the constraint Q — Qi. These 
constrained stationary densities were then improved upon iteratively. The details of p 1 / 2 jr, the counterpart of ip of 
Eq. (0), are obtained more precisely when the r 2 behaviour near r = and the harmonic oscillator asymptotics at 
infinity are factored out in Eq.(H). In numerical work, we use a function a{r,r) such that p(r,r) = r 2 e~ r e 2a ^ r ' T \ 
and properly express the term o^f — d r ((d T f) 2 /p) - see Appendix. 

We also performed the minimization of the functional (|9|) treating p(rj,Qi) as independent variables. It turns out 
that p(r, t) thus obtained do not fulfil Eq.(||) accurately, but the accuracy in S is better than 0.1%. 

The numerical results have been obtained with physical data on 7 Li adopted after the most accurate treatment 
up to date j7j. These values give K = —5.74 x 10~ 3 x N, and we obtain the critical value K c between -7.2249 and 
-7.2255 ( N c between 1258.7 and 1258.8). 

The potential energy V(Q) from the bounce solutions (Fig.l), nearly identical with V[p s {Q)] for constrained sta- 
tionary p s (r, Qi), is very flat between Q(0) and Q(t p /2) for N N c . For smaller N, it becomes quite peaked around 
the summit, and its fall on the side of small Q becomes very steep. All the obtained bounce solutions result from 
the small adjustment of the initial densities. For larger N c — N (and increasing Q 2 and Q terms, cf Eq.(A4)) this 



adjustment becomes gradually more difficult. We could not obtain the exact solution of Eq.rtq) (or even a constrained 
stationary p s (r, Qi) for small Qi) for TV < 1200. 

The mass parameters B(Q) from various instanton solutions (Fig. 2) are nearly identical which shows that there 
exists a universal inertia for the radial collective motion of BEC close to instability. This may be understood as a 
consequence of a nearly static character of solutions f(r,Q) in the limited range of N(K) values of interest: / depends 
weakly on N close to critical N c . 

For a gaussian density with a variable width b(r), p(r,r) — n~ 3 / 2 b~ 3 r 2 exp(— r 2 /b 2 ), the current j — b(r/b)p 
(Eq.(|)), and the mass parameter B(b) = Air / b 2 J °° r 2 pdr = 3/2 (cf P,|i"4|). Hence, as for the gaussian density 
Q = 3b 2 /2, one has AQB{Q) = B(^/Q) = 1, as used in Ref. 0. As seen in Fig. 2, this is a fair assumption near 
the metastable minima, much worse though for smaller Q around the barrier's summit. The error in S due to the 
gaussian value of B amounts to 3-4% in the cases studied. 

The bounce "amplitude", ^fpjr (Fig. 3), differs from the constrained stationary values ^fpljr mostly near r = and 
for small Q, by up to 1%. Decay exponents calculated with the initial densities p s {r,Qi) are only up to 0.3% larger 
than the exact ones. 

Next we turn to collective radial oscillations of BEC and their decay rates. Quantum tunneling from the excited 
state with energy hco n may be treated by solving Eq.(^) as before, except that now the period t p must be finite and 
N£ replaced by N£ n = N£ + fuv n . One has to fix the boundary values p(r, t p /2) and the energies Tiw n . 

The energies hu) n of the excited quantum states follow from the quantization condition S = 2mr p5| , where 

S = 2N Jq^ /2) dQ v / 2B(Q){£ n - V(Q)), cf Eq.(|), with the period t p depending on %u n . Essentially, the imaginary- 
time boundary value p(r,T p /2) should match the real-time oscillatory solution p(r,t p /2) at the turning point Q at 
which V(Q) =£ n . 

We could not solve the real-time version of Eq.(||) as accurately as for instantons. However, since the action S is 
insensitive to the details of solution, the lowest quantized hu) n are quite accurate. Their comparison to frequencies of 
the small amplitude oscillations, obtained from the GPE linearized around ipo (so called Bogolyubov spectrum, see 
Jl6|]), is interesting. The first radial excitation, tkoi, is nearly equal, but little smaller than, the lowest Bogolyubov 
mode corresponding to the small amplitude limit, except that it may not exist for a too shallow V(Q), like for 
N = 1255. The higher modes, TiLo n , n — 2,3,..., come out roughly at the multiples of Tlloi (Table 1), and thus 
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represent nearly harmonic spectrum of collective radial oscillations. Its n — 2,3,4 states are lower than the second 
Bogolyubov mode, lying slightly above AhajQ. 

Using the quantized oscillation energies we have found periodic instantons and calculated decay exponents for some 
low collective radial excitations. Such instanton with a period r p is quite similar to the so-called "thermon" which 
describes thermal decay from the metastable ground state at the temperature kT — h/r p [^|. This similarity is easily 
understood: Thermal decay goes via thermal excitation of quantum states above the metastable ground state and 
successive tunneling from them through the smaller barrier. At higher temperatures, details of the excitation spectrum 
and tunneling become irrelevant and the thermal tunneling rate exp(— (Eb — E)/(hujo)), with Eb = N£b the barrier 
energy, replaces a thermal mixture of rates exp(— S n ) from a few lowest excited states at small temperature. The 
critical temperature T c at which this happens is usually determined as the one at which T t hermai(T c ) = r t ^ ermon (r c ). 

Numerical results are collected in Table 1. Those in columns 2-4 depend only on K and relate to all BEC-s with 
attractive interaction. Decay rates T from the metastable 7 Li BEC were calculated for ujq — 908.41 s _1 g], assuming 
prefactor (T/u;o)e s = (uji/luo)(15S/2'k) 1 / 2 with U\ from the small amplitude limit. The exact prefactor is difficult 
to calculate, but it must be of the same order as in Table 1. Obtained values of T are very close to those of Ref. 0, 
where the effect of too small B-s is fortuitously reduced by half by the interpolation overestimating V(Q). 

The decrease of the decay exponents S n with the mode number n and their comparison to exponents S for the 
ground state are seen in Table 1, col. 5. For example, quantum decay of the n = 2 radial mode in BEC with 
TV = 1245 is nearly as quick as for the ground state of the N — 1255 BEC. The crossover temperature T c at 
which the thermal decay begins to dominate over the quantum tunneling is approximated, neglecting prefactors, as 
kT c /(huJo) — (Eb — E)/S, with Eb = N£b the barrier energy, E — N£ for the metastable state, and E = N£ n 
for the excited state. For the latter, the meaning of critical temperature is extended in analogy with that for the 
ground state: Suppose the condensate in the oscillation state n and ask at which temperature T c (n) the thermal rate 
exponent equals the instanton action S n . These ratios (Table 1, column 3) show that for presently measurable T and 
wo = 908.41 s _1 , corresponding to T = 6.94 nK, T c w InK, both for ground state and collective radial oscillations. 
We notice, that for all finite- t p instantons describing decay out of the n— th collective radial excitation, Tl/t p < T c (n) 
so, at their own thermon "temperature", they dominate over the thermal decay 

The quantity S/N (Table 1, column 4) shows (N c - N)€ behaviour close to N c , with £ = 5/4 ffift. For larger 
N c — N the effective exponent slightly decreases to £ « 1.2. Consider now two different attractive BEC-s with critical 
particle numbers N c and N' c . The same decay exponent S is obtained for such (not too large) N c — N and N' c — N' 
which satisfy the relation N' c — N' = (N c - N)(N^/N C ) 1 - 1 ^, with 1 - 1/f « 1/5. Thus, e.g. we obtain S = 9.44 for 
N[.-N' w 12 in BEC with Nf. = 6294, and for N' c - N' = 6 - 7 in BEC with N' c = 251.76, cf Table 1. 

In summary, the equations for the condensate density, describing both the real- and imaginary-time dynamics of 
spherical BEC, were formulated and the exact instanton solutions were found numerically, also for collective radial 
excitations. The determined mass parameter (Fig. 2) deviates from the gaussian ansatz, but the calculated decay 
exponents for the metastable states agree well with Ref. 0. It follows from Fig. 2 that the exact mass parametr 
may be more relevant to the behind-the-barrier collapse phase, i.e for smaller Q. The quantized energies of collective 
finite amplitude radial vibrations form nearly harmonic (slightly compressed) spectrum with uj n sa nwi, where u>i is 
slightly lower than the lowest Bogolyubov mode. The n = 2, 3, 4 collective oscillation states lie lower than the second 
Bogolyubov mode. Any excitation (thermal or otherwise, e.g. by modulation of the trapping oscillator frequancy) of 
these collective vibrations must lead to a faster decay of the condensate, as the Table 1 shows. 

If quantum tunneling is not to be overshadowed by thermal decay, the experiments should proceed at low T and/or 
large loq. Since theoretical results are more certain away from N c , where the exponent S dominates decay, one should 
probe a range of moderate S, giving observable, but not too large T (perhaps, by discarding prompt collapses). The 
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corresponding range of N c — N values depends on N c as N c ■ 



APPENDIX: NUMERICAL METHODS 



We have p(s, t) 



= se 



s e 2a(s,r) s _ r "2. _ ^phe stationary GPE leads to the equation for a(s) 



with (3 — e — 3/2. For large s, a 




(Al) 



(A2) 



A solution regular at s = must fulfil 



4 



§(0) = |(3/2 -e + Ke^). 



(A3) 



These boundary conditi ons suggest a method of soluti on: For a given K we assume som e e and C and, starting from 
the asymptotic values (A2) at large s, integrate Eq.(Al) to s = 0. We check Eq.(A3) and the normalization and 
correct e and C until we fulfil both. 

By using a new variable v = ^/ (rp) and factoring out 2rp, we transform Eq. (||) to a form 



9a 



2dQ dQ 



+ v 2 [l + s{2 



da 



1)]) 



dR 
7to' 



(A4) 



suitable for both instantons and oscillati ons if Q 2 = 2(V(Q) — £„)/B(Q) is understood. R[a] stands for the l.h.s. 
of E q. ( Al| ) . Let us call the l.h.s. of Eq.(A4) F. If F = (no time dependence) we recover stationary solutions of 
Eq.(AlJ7~For instantons, we solve Eq.(A4) iteratively. Having a set of a(s,Qi) we calculate Fi = F(s,Qi). For each 
Fi we solve (A4) as the ordinary differential equation in s to obtain new a(s,Qi). The method, as for the stationary 
case, is to adjust the asympotic form (A2) to the proper regularity condition at s — 0. The new and old Fi-s are 
combined to provide initial Fi-s for the next iteration. With a careful modification of Fj-s this iteration leads to the 
self-consistency, i.e. Fi(old) = Fi(new). The initial densities p s (s,Qi) are obtained using the cons trai ned imaginary- 
time step Hartree procedure. We use the Runge-Kutta-Merson procedure for integration of Eqs.(A4). Energies are 
calculated using a mesh of 128 equidistant points, r/do = — 4.5, and the cubic spline interpolation for derivatives. 
We have checked that doubling the mesh density does not change results in any appreciable manner. 
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Table 1 - Energies, crossover temperatures, decay exponents and rates for metastable and radially excited states of 
the 7 Li BEC. Results from bounce solutions (* from the functional minimization). 



N 




kT c 


x 10 3 

N 


S -S-e 6 ' 


L J 




1ZOO 




n 1 ^o 

U. loU 


9 P.7AA 
Z.u 1 44 


O.OOO 4.OU0 


i "zr 1 n 

l.oD -1U 


•2 


1 o^n 

IZOU 




n 1 


7 ^90 


O AA(\ A 7Kfi 
y .44U 4. / DO 


Q /I Q 10" 
o.4o -1U 


-I 


n=l 


0.967 


0.188 


2.5244 


3.156 










n 1 qo 
U. lot) 


1 Q H07 


16.306 6.87 


5.17 10" 


-4 


11 — 1 


l.UoO 


n 9ns 

U.ZUo 


7 RQ'i/l 

/ .oyo4 


9.578 






, n 9 

n — z 


9 119 
Z. 1 1Z 


n 91 7 

U.Zl / 


o.OOou 


4.437 






1 940 




n 907 


1 Q 089 


23.662 8.80 


4.23 10" 


-Y 


n=l 


1.169 


0.224 


13.436 


16.661 






n=2 


2.292 


0.234 


9.0346 


11.203 






n=3 


3.373 


0.241 


5.1464 


6.382 






1230 




0.238 


31.980 


39.335 12.31 9.24 -10" 


14 


n=l 


1.268 


0.252 


26.084 


32.083 






n=2 


2.518 


0.260 


21.347 


26.257 






n=3 


3.748 


0.267 


17.076 


21.004 






1200* 




0.305 


75.40 


90.490 21.12 9.63 -10" 


36 



Figure captions 

Fig.l Potential energy £(Q) = NV(Q) of BEC (in Tiwo) for various N < N c . (Q in units d§) 

Fig. 2 Mass parameters J5( v / Q) = AQB{Q) from various instanton solutions, overlayed in one picture. For gaussians, 
4QB(Q)=1. 

Fig. 3 Bounce penetrates the barrier practically in a finite r (in units to^ 1 ). The metastable density p(r, ±oo) is 
nearly equal to p(r, r = 3.43) shown. 
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